In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.optimize import curve_fit

Models


In [100]:
def gauss(x, a, center, sigma, offset):
    return a * np.exp(-(x-center)**2/sigma**2) + offset
def twogauss(x, a1, center1, sigma1, offset1, a2, center2, sigma2, offset2):
    return gauss(x, a1, center1, sigma1, offset1) + gauss(x, a2, center2, sigma2, offset2)
def gausslinear(x, a, center, sigma, offset, slope, intercept):
    return gauss(x, a, center, sigma, offset) + slope * x + intercept
def twogausslinear(x, a, center, sigma, offset, a1, center1, sigma1, offset1, slope, intercept):
    return gausslinear(x, a, center, sigma, offset, slope, intercept) + gauss(x, a1, center1, sigma1, offset1)

Single Gaussian Fit


In [34]:
number_of_events = 1000
bin_range = [-10, 10]
bin_width = 0.3

In [35]:
bins = np.arange(bin_range[0], bin_range[1], bin_width)
bin_mids = np.arange(bin_range[0], bin_range[1], bin_width) + bin_width/2.

In [69]:
line_strength = [0.4, 0.6]
line = np.random.randn(number_of_events*line_strength[0])
background = np.random.uniform(bin_range[0], bin_range[1], size=number_of_events*line_strength[1])
photons = np.append(line, background)

In [70]:
hist, bin_edges = np.histogram(photons, bins = bins)

In [71]:
ydata = hist
xdata = bin_mids[:-1]
plt.bar(bin_edges[:-1], hist, width = bin_width);
plt.xlim(min(bin_edges), max(bin_edges));
plt.plot(bin_mids[:-1], hist)


Out[71]:
[<matplotlib.lines.Line2D at 0x108103cd0>]

Fit


In [72]:
ydata = hist
xdata = bin_mids[:-1]
fit_guess = [1, 0, 1, 0, 1, 1]
popt, pcov = curve_fit(gausslinear, xdata, ydata, sigma=np.sqrt(ydata)+1, p0=fit_guess)

In [73]:
popt


Out[73]:
array([  5.07269116e+01,   2.44597480e-02,  -1.42691635e+00,
         4.14966363e+02,  -8.54953992e-02,  -4.06930553e+02])

In [77]:
plt.bar(bin_edges[:-1], hist, width = bin_width);
plt.xlim(min(bin_edges), max(bin_edges));
plt.ylim(0)
#plt.plot(bin_mids[:-1], hist)
plt.plot(xdata, gausslinear(xdata, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5]), label='Fit')
plt.legend()


Out[77]:
<matplotlib.legend.Legend at 0x10b832190>

In [75]:
perr = np.sqrt(np.diag(pcov))

In [76]:
perr


Out[76]:
array([  2.87543297e+00,   5.30219604e-02,   7.19252050e-02,
         2.36600951e+07,   5.81707460e-02,   2.36600951e+07])

In [101]:
number_of_trials = 10000
centers = np.zeros(number_of_trials)
centers_errors = np.zeros(number_of_trials)

def random_gauss_fits(number_of_trials):
    
    line_strength = [0.4, 0.6]

    for i in np.arange(0, number_of_trials):
        line = np.random.randn(number_of_events*line_strength[0])
        background = np.random.uniform(bin_range[0], bin_range[1], size=number_of_events*line_strength[1])
        photons = np.append(line, background)
        hist, bin_edges = np.histogram(photons, bins = bins)
        fit_guess = [1, 0, 1, 0, 1, 1]
        popt, pcov = curve_fit(gausslinear, bin_mids[:-1], hist, sigma=np.sqrt(ydata)+1, p0=fit_guess)
        centers[i] = popt[1]
    return centers

In [102]:
centers = random_gauss_fits(number_of_trials)

Distribution of Fit Centers


In [110]:
plt.hist(centers, bins=20, label='$\sigma$ = ' + "%0.3f" % centers.std(), range=[-0.5,0.5])
plt.xlabel("Fit center")
plt.xlim(-0.5, 0.5)
plt.legend()


Out[110]:
<matplotlib.legend.Legend at 0x10d7b0d50>

Two Gaussians


In [130]:
number_of_events = 1000
line_strength = np.array([0.4, 0.4, 0.2])
line_positions = [0, 10]
bin_range = [-10, 20]
bin_width = 0.3

In [131]:
bins = np.arange(bin_range[0], bin_range[1], bin_width)
bin_mids = np.arange(bin_range[0], bin_range[1], bin_width) + bin_width/2.

In [132]:
line1 = np.random.randn(number_of_events*line_strength[0]) + line_positions[0]
line2 = np.random.randn(number_of_events*line_strength[1]) + line_positions[1]
background = np.random.uniform(bin_range[0], bin_range[1], size=number_of_events*line_strength[2])
photons = np.append(line1, np.append(line2, background))

In [133]:
hist, bin_edges = np.histogram(photons, bins = bins)

In [134]:
plt.bar(bin_edges[:-1], hist, width = bin_width);
plt.xlim(min(bin_edges), max(bin_edges));
plt.plot(bin_mids[:-1], hist)


Out[134]:
[<matplotlib.lines.Line2D at 0x10db8b250>]

In [10]:
popt, pcov = curve_fit(twogausslinear, xdata, ydata, sigma=np.sqrt(ydata)+1, p0=[1, 0, 1, 0, 1, 10, 1, 0])

In [11]:
popt


Out[11]:
array([  9.50964076e+01,  -2.44937032e-02,   1.38935303e+00,
         7.99067997e+01,   2.59030560e+01,   1.00315649e+01,
         1.24359943e+00,  -7.99520670e+01])

In [12]:
pcov


Out[12]:
array([[  6.11554534e+00,  -1.29398291e-03,  -3.47002223e-02,
          1.78619449e+04,   1.63665243e-03,   7.46230510e-06,
         -8.12015335e-04,  -1.78619393e+04],
       [ -1.29398294e-03,   5.10982021e-04,   2.03100396e-05,
          1.40794303e+03,   4.70732715e-06,   2.57796229e-08,
         -2.19259263e-06,  -1.40794302e+03],
       [ -3.47002223e-02,   2.03100384e-05,   6.02941493e-04,
          3.11669352e+03,  -7.98020166e-05,  -3.51174601e-07,
          4.00132563e-05,  -3.11669382e+03],
       [  1.78626616e+04,   1.40793974e+03,   3.11665515e+03,
         -3.78035305e+13,  -8.80817424e+03,  -1.60122181e+02,
          3.98058125e+02,   3.78035307e+13],
       [  1.63665261e-03,   4.70732710e-06,  -7.98020240e-05,
         -8.80819344e+03,   2.25888267e+00,   6.86695748e-03,
         -4.48027080e-02,   8.80819480e+03],
       [  7.46230819e-06,   2.57796131e-08,  -3.51174754e-07,
         -1.60122287e+02,   6.86695748e-03,   2.23965485e-03,
         -3.95895585e-04,   1.60122294e+02],
       [ -8.12015343e-04,  -2.19259267e-06,   4.00132565e-05,
          3.98060051e+02,  -4.48027080e-02,  -3.95895585e-04,
          2.76672551e-03,  -3.98060704e+02],
       [ -1.78626560e+04,  -1.40793971e+03,  -3.11665542e+03,
          3.78035307e+13,   8.80817560e+03,   1.60122188e+02,
         -3.98058777e+02,  -3.78035310e+13]])

In [13]:
plt.bar(bin_edges[:-1], hist, width = bin_width);
plt.xlim(min(bin_edges), max(bin_edges));
plt.ylim(0)
plt.plot(bin_mids[:-1], hist)
plt.plot(xdata, func(xdata, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6], popt[7]), label='Fit')
plt.legend()


Out[13]:
<matplotlib.legend.Legend at 0x10780f490>

In [17]:
pcov


Out[17]:
array([[  6.11554534e+00,  -1.29398291e-03,  -3.47002223e-02,
          1.78619449e+04,   1.63665243e-03,   7.46230510e-06,
         -8.12015335e-04,  -1.78619393e+04],
       [ -1.29398294e-03,   5.10982021e-04,   2.03100396e-05,
          1.40794303e+03,   4.70732715e-06,   2.57796229e-08,
         -2.19259263e-06,  -1.40794302e+03],
       [ -3.47002223e-02,   2.03100384e-05,   6.02941493e-04,
          3.11669352e+03,  -7.98020166e-05,  -3.51174601e-07,
          4.00132563e-05,  -3.11669382e+03],
       [  1.78626616e+04,   1.40793974e+03,   3.11665515e+03,
         -3.78035305e+13,  -8.80817424e+03,  -1.60122181e+02,
          3.98058125e+02,   3.78035307e+13],
       [  1.63665261e-03,   4.70732710e-06,  -7.98020240e-05,
         -8.80819344e+03,   2.25888267e+00,   6.86695748e-03,
         -4.48027080e-02,   8.80819480e+03],
       [  7.46230819e-06,   2.57796131e-08,  -3.51174754e-07,
         -1.60122287e+02,   6.86695748e-03,   2.23965485e-03,
         -3.95895585e-04,   1.60122294e+02],
       [ -8.12015343e-04,  -2.19259267e-06,   4.00132565e-05,
          3.98060051e+02,  -4.48027080e-02,  -3.95895585e-04,
          2.76672551e-03,  -3.98060704e+02],
       [ -1.78626560e+04,  -1.40793971e+03,  -3.11665542e+03,
          3.78035307e+13,   8.80817560e+03,   1.60122188e+02,
         -3.98058777e+02,  -3.78035310e+13]])

In [14]:
perr = np.sqrt(np.diag(pcov))


/Users/schriste/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:1: RuntimeWarning: invalid value encountered in sqrt
  if __name__ == '__main__':

In [15]:
perr


Out[15]:
array([ 2.47296287,  0.02260491,  0.02455487,         nan,  1.50295797,
        0.04732499,  0.05259967,         nan])

In [16]:
x = 
plt.plot()


  File "<ipython-input-16-ebb8e0862158>", line 1
    x =
        ^
SyntaxError: invalid syntax

In [ ]: